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I.  INTRODUCTION 


The  stability  of  ship  towing  is  the  primary  concern  in  both  naval  architecture  and 
maneuvering  control.  The  main  problem  consists  of  two  individual  tasks,  first  building  a 
mathematical  model  of  the  two  towed  ships  and  second  analyzing  the  directional  stability 
of  the  combined  system. 

Traditionally,  ship  towing  in  the  open  ocean  has  been  performed  under  relatively 
large  separation  distances  between  the  two  vessels;  i.e.,  large  towlines.  Only  in  the  case 
of  barge  towing  in  rivers  and  canals,  has  close  proximity  towing  been  applied.  Over  the 
last  few  years,  however,  interest  on  ocean  close  proximity  towing  has  been  resurfaced. 
The  Office  of  Naval  Research  is  interested  in  close  proximity  towing  for  its  Advanced 
Hull  forms  Program,  in  particular  small  waterplane  area  twin  hull  (SWATH)  ships  and 
their  variations  (such  as  the  SLICE  hull).  Studies  of  the  SEA  LANCE  project  at  the  Naval 
Postgraduate  School  (Total  Ship  Systems  Engineering  Program)  have  also  shown  several 
benefits  associated  with  close  proximity  towing.  One  of  the  main  benefits  of  close 
proximity  towing  in  military  applications  is  the  ability  to  separate  a  combatant  ship  from 
a  main  part  of  its  payload. 

Towing  is  not  without  its  problems,  however.  Directional  stability  under  tow  and 
seakeeping  are  two  main  areas  of  concern.  This  thesis  will  concentrate  on  the  problem  of 
directional  stability.  Previous  studies  on  directional  stability  of  ship  towing  were 
performed  by  Abkowitz  in  1964  who  developed  the  characteristic  equation  for  single 
body  towing,  and  by  Bemitsas  and  others  who  developed  the  criteria  for  stability  of 
Abkowitz’s  4//!  order  characteristic  equation.  In  these  studies  the  stability  criteria  was 
based  solely  on  the  trailing  ship,  which  means  that  the  dynamics  of  the  leading  ship  were 
neglected.  In  fact,  the  leading  ship  was  approximated  by  a  point  mass  moving  on  a 
straight  line  under  constant  forward  speed.  Our  concern  is  that  although  this  may  be  valid 
for  long  towlines,  it  may  not  be  a  valid  approximation  for  close  proximity  towing,  and  the 
existing  stability  criteria  may  be  inadequate  which  is  the  reason  and  the  justification  of 
this  thesis. 
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In  this  study  as  a  leading  ship  we  will  use  the  SLICE  vessel,  which  is  105  ft,  and 
180  tons,  and  as  a  trailing  ship  the  SWATH  ship  Kaimalino  which  is  89  ft,  and  217  tons. 

The  Slice  is  a  high  speed  variant  of  the  Swath  technology.  It  has  4  underwater 
hulls  instead  of  two.  Attached  to  each  hull  is  a  strut  that  extends  up  to  support  the  main 
body.  Figure  1  shows  the  profde  view  of  the  Slice.  [2] 


The  SSP  Kaimalino  was  the  world’s  first  high  performance  open  ocean  Swath 
ship.  It  consists  of  two  parallel  torpedo-like  hulls.  Attached  to  the  hulls  are  two 
streamlined  struts.  The  struts  extend  above  the  water  surface  and  support  the  main  body. 
The  Kaimalino  also  has  stabilizing  fins  attached  near  the  aft  end  of  each  hull.  Figure  2 
shows  a  profile  of  the  SSP  Kaimalino.  [2] 


Figure  2.  The  profile  of  the  SSP  KAIMAFINO  [From:  [2]] 
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The  general  approach  adopted  in  this  thesis  is  as  follows: 

1.  First  we  will  use  the  maneuvering  equations  of  motion  in  the  horizontal  plane  for 
both  leading  and  trailing  ship. 

2.  Coupling  between  the  two  ships  will  be  provided  through  the  towline. 
Hydrodynamic  coupling  arising  from  radiation  and  diffraction  effects  will  be 
neglected.  This  can  be  easily  incorporated  into  the  analysis  once  the  effect  of  such 
hydrodynamic  coupling  on  the  hydrodynamic  coefficients  is  established. 

3.  The  response  of  the  coupled  system  will  be  analyzed  by  both  simulation  and  an 
eigenvalue  analysis.  It  is  hoped  that  this  analysis  will  provide  regions  of 
directional  stability,  and  can  therefore  lead  to  a  design  methodology  for  rational 
towing  system  parameter  selection  based  on  their  stability  properties. 

The  hydrodynamic  coefficients  of  both  ships  used  in  this  study  are  provided  from 
reference  [2],  It  should  be  emphasized  that  our  procedure  will  be  general  enough  so  that  it 
can  accommodate  different  or  additional  hydrodynamic  coefficients. 
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II.  PROBLEM  FORMULATION 


A.  EQUATIONS  OF  MOTION 

The  surge,  sway  and  yaw  maneuvering  equations  of  motion  are  shown  below. 
Throughout  this  thesis  subscript  1  is  for  the  towing  (leading)  ship  and  subscript  2  is  for 


the  towed  (trailing)  ship. 

(m2  -X.2)u2  =  -mv2r2  -R2+T cos(^2  +  y)  (1) 

(m2  ~  Ki  )v2  -  y,-2>\  =  ~m2r2u2  +  Yv2v2  +  Yr2r2  -  T  sin (y/2  +  y)  (2) 

(1 Z2  ~  ^ rl)^2  ~-^v2Y2  ~  ^ v2V2  T  ^ r2r2  ~  Tx  p2SVd(\j/2  +/)  (3) 

(m,  -  Yn  )vj  -  Ynrx  =  -mxrxux  +  YvXvx  +  YrXrx  +  T  sin(^  +  y)  (4) 

(4i  -  Nn  )?'i  -  =  NaY  +  NX\  ~  TxPi  sin(^i  +  7)  (5) 

and 

Vi  =  r\  (6) 

V  2  =  r2  (7) 


The  detailed  explanation  of  the  derivation  of  the  equations  of  the  motion,  and  the 
assumptions  used  can  be  found  in  the  reference  [1]. 

In  these  equations  u  is  the  surge  velocity,  v  is  sway  velocity,  r  is  the  yaw  velocity, 
R  is  the  resistance  of  the  vessel  moving  through  body  of  water,  and  T  is  the  tension  in  the 
connection  (rope,  cable,  or  some  other  mechanism)  between  the  two  towing  ships. 

The  kinematic  relations  are  as  shown  below; 

yx  =ux  sin  y/x  +  vx  cos  y/x  (8) 

y2  =  u2  sin  y/2  +  v2  cos  y/2  (9) 

If  we  define 
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y  —  y i  ~y\ 


(10) 


then 

y  =  y2-y  \  (H) 

These  geometric  relations  are  explained  in  Figure  3. 


Substituting  equation  (8)  and  equation  (9)  into  equation  (10),  we  get 
y  =  u2  sin  y/2  +  v2  cos  y/2  -ul  sin  y/x  -  vi  cos  y/] 

From  the  geometry  of  the  figure  we  can  see  that; 
y2  +  x  2  sin  ^2  -  (yl  -  x  x  sin  y/x ) 

nn-i  /IX  -  r  r 


which  can  be  rewritten  as 


smy  ■ 


v  1  ( 

-  +  -\xp2  sm y/2  +x 


/  / 


p  i 


(12) 


(13) 


(14) 
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B.  LINEARIZATION 

In  order  to  analyze  the  stability  of  our  system,  we  should  first  linearize  the 
system.  During  linearization  we  assumed  that  the  velocities  of  the  both  vessels  are 
constant  and  identical  which  in  non-dimensional  terms  means 

U\  =  M2  =  1  (15) 

Another  assumption  is  that  we  are  dealing  with  small  heading  angles,  which  gives 
us 

sin  ¥\  ,i=y/i  ,2  (16) 

and 

cos^j  2  =  1  (17) 

Since  we  are  dealing  with  constant  surge  velocity,  equation  (1)  is  dropped,  and  we 
are  left  with  the  equations  (2),  (3),  (4),  and  (5). 

When  we  substitute  equations  (15),  (16),  (17)  into  equations  (2),  (3),  (4),  (5),  (12) 


and  (13)  we  come  up  with 

(in 2  -  Yv2  )v2  -  Yrli'2  =  ~m2r2  +  Yv2V2  +  Y r2>\  ~  T (V 2+7)  (18) 

(J Z2  —  N j.2)Y2  —  ^v2+2  ~Nv2V2  +  ^ rir2  ~  YX  p2(l// 2  +  7)  (19) 

K  -  Yn )vx  -  Ynrx  =  -m,r,  +  Fvlv,  +  YrXrx  +  T(i//X  +  y)  (20) 

(4i  -  Nn  )fi  -  Nnvx  =  Nvlvx  +  NrXrx  -  TxpX  (if/,  +  y)  (21) 

y  =  ¥ 2  +  v2  ~V\~ V1  (22) 

r=|  +  y(^2^ 2+xpM)  (23) 

The  above  equations  can  be  rewritten  into  a  standard  matrix  form  as  follows: 

^2  —  ^2vvV2  +  ^2vrr2  +  ^2v  (^2  +  7 )  (24) 
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2)  4nA2  A2ir^2  ^  4.  (l4  y)  (25) 

>)  =  4,'i  +  4,,'i  +  B„  (v,  +  r)  (26) 

'i  =  4,,1’,  +  Ai„n  +  B„  (y,  +  r)  (27) 

After  the  necessary  mathematical  steps  we  can  get  the  coefficients  as  follows 

4„  =  [(2.2  -A4K2  +(AT,!lrJM('»2 -y»ll,2 -N„)-N„Yn\  (28) 

4„  =[()■■.: -»nl(V-,-v(j-(v(A;..)]-[(/«, Af,r„]  (29) 
4,  =[-(2,2 -Nt2)T-[YnTxp2M(m2 -Yn\lzl -Nn)-NnYn]  (30) 

4„,  =[fcA„)+ A,2(m2  -  f  ,  )p[(m2  -  -  K / - ,  -  V).  )-NtlYn]  (31) 

/!.„=  [(!',, -m.I.V,,  +Nrl(m1-Yn)\+[(m1-Ynilzl-Nn)-NnYn]  (32) 

B2r=\-TNll-Txp2(m1-Yn)\+[(m1-Ynll,1-Nn)-NnYn\  (33) 

4„  =[(2,.  -A2„)i;,  +(M„4)]-[(™,  -n,X/„  -2Vrl)-Ma7j  (34) 

4..  =  [(4  -m,)(2.,  -4, )+((v,4, ,)]+[('»,  -n,X4,  -JV„)-JV„r„]  (35) 

K  =  [(l.2-Ntt)T-[YnTxrlM(mt-Y„llzl-Nn)-N„Yn]  (36) 

4„  =[fc,(Vn)+(V„('»i-4,)]-[('»,  -4X2., -JVj-JV.,7,,]  (37) 

4 -'X.K,  +  2Vr, (m.  -4,)MK -7„X2.1  -Ar„)-Ar„y„]  (38) 

Btr=[TN„  +Txpl(m, -iyH(<n, -Y^I  zl -Nn)-N„Yn\  (39) 

Now  we  can  write  down  the  complete  system  in  a  7x7  matrix  form  as 
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^2 

v2 

r2 

A 

A 

fi 

=  U] 

>\ 

y 

y 

¥2 

¥2 

¥x_ 

where 


U] 


^2w  ^2  w  0  0 

^2  rv  Alrr  0  0 

0  0  Alvv  Alvr 

0  0  Alrv  Alrr 

10-10 
0  10  0 
0  0  0  1 


B 


2v 


B , 


B2v  +  ■ 


B-x 


2vJvp2 


B2  r  + 


1 

B2rXp2 


/ 

B\rXp2 


l 

K  Kx 

i 

K 

i 

o 
o 
0 


pi 


l 


1 

0 

0 


B2vX  p\ 


B2rXp\ 


B  iv  + 


B.x 


Blr+- 
-1 
0 
0 


hWpl 

/ 


I 


(40) 


(41) 


Following  linearization  of  the  system,  the  next  step  is  to  find  the  eigenvalues  of 
the  A  matrix.  These  will  determine  whether  the  system  is  stable  or  not.  The  matlab 
program  shown  in  Appendix  D  was  written  and  it  was  used  to  find  the  eigenvalues  of  the 
system. 

C.  A  ZERO  EIGENVALUE 

Several  different  values  for  tension  and  for  length  of  the  connection  between  the 
two  ships  were  tried  using  the  program  shown  in  Appendix  D.  In  every  instant,  a  zero 
eigenvalue  was  found.  Therefore,  it  seemed  impossible  for  the  system  to  reach  a  stable 
condition.  Simulations  also  showed  the  same  effect.  In  order  to  test  whether  the  existence 
of  a  zero  eigenvalue  was  just  a  coincidence  or  an  inherent  property  of  the  system,  the 
matlab  program  shown  in  Appendix  E  is  used.  This  performed  a  symbolic  manipulation 
utilizing  the  Matlab  symbolic  manipulation  toolbox.  It  was  found  that  regardless  of  the 
numerical  values  of  the  geometric  properties  and  the  hydrodynamic  coefficients  of  the 
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system,  the  characteristic  equation  of  this  system  has  always  a  zero  constant,  which 
proved  the  existence  of  a  zero  eigenvalue.  Physically,  this  zero  eigenvalue  can  be 
explained  by  the  lack  of  any  directional  stabilizing  effects  on  the  leading  ship.  It  can  be 
easily  verified,  by  observing  the  signs  of  the  towline  restoring  forces  and  moments  on  the 
two  ships  that  the  towing  force  has  a  stabilizing  effect  on  the  trailing  ship  but  a 
destabilizing  effect  on  the  leading  ship.  Since  ships  in  the  horizontal  plane  cannot  exhibit 
directional  stability,  it  is  necessary  in  order  to  continue  with  the  analysis  to  introduce 
some  kind  of  directional  control  on  the  leading  ship.  This  is  analyzed  in  the  following 
section. 

D.  RUDDER  CONTROL 

It  was  shown  in  the  previous  section  that  some  kind  of  directional  control  on  the 
leading  ship  is  necessary  to  proceed  with  the  towing  analysis.  We  do  this  my  means  of  a 
rudder  control  law.  The  equations  of  motion  (4),  (5)  are  changed  because  of  inclusion  of 


the  rudder,  but  equations  (2)  and  (3)  are  same. 

(*n2  -  Yvi )v2  -  Yf2l%2  =  -m2r2u2  +  Yv2v2  +  Yr2r2  -  T  sin (y/2  +  y)  (2) 

(I Z2  ~  ^  i-2  _  ^ vjY 2  N. v2V2  r 2^2  ~  Yx  p2  Y )  O) 

(m,  -  Yvt  )vj  -  Ynrx  =  -mxi\ux  +  Yvlvx  +  Yrlrx  +  T  sin(^,  +y)  +  YsS  (42) 

(4i  -  Nn  )fi  -  NvA  =  Nvivi  +  Nriri  -  Txpi  sin(^i  +  r)  +  NsS  (43) 

where 

8  =  ~kVWx  -  KYx  -  krr\  ~  ky}’\  (44) 


The  rudder  control  law  shown  in  equation  (44)  is  a  typical  full  state  feedback 
control  law.  This  was  chosen  as  a  representative  control  law  and  it  by  no  means  limits  the 
applicability  of  the  results  that  follow.  The  procedure  would  be  identical  even  if  a 
different  rudder  control  law  were  employed.  The  feedback  gains  were  calculated  by  using 


10 


the  standard  pole  placement  scripts  in  Matlab’s  control  system  toolbox.  In  order  to 
present  the  results  in  a  compact  form,  we  chose  the  control  time  constant  as  the 
representative  selection  criterion  for  the  gains.  Different  choices  do  not  alter  the 
qualitative  features  of  out  results.  The  control  time  constant  is  defined  as  a  non- 
dimensional  response  time  for  the  system.  It  is  measured  in  dimensionless  seconds,  with 
one  dimensionless  second  being  the  time  that  it  takes  for  a  ship  to  travel  one  ship  length. 
Typically,  a  system  reaches  steady  state  within  three  time  constants.  Thus,  a  high  time 
constant  signifies  a  rather  slow  leading  ship  control  law,  while  a  small  time  constant 
shows  a  fast  (responsive)  control.  The  closed  loop  poles  are  simply  set  as  the  negative 
inverse  of  the  time  constant. 

E.  SAMPLE  SIMULATION  RESULTS 

We  want  to  give  one  example  to  each  stable  and  instable  condition  to  clarify  some 
of  the  concepts  mentioned  above.  These  results  were  obtained  using  the  matlab  code 
shown  in  Appendix  H. 

1.  Stable  Condition 

The  figures  below  are  obtained  for  the  following  parameters; 

Tension=0.005 
Length=1.2 
Time  Constant=l 

The  results  show  a  slow  convergence  to  the  nominal  equilibrium  state;  i.e., 
straight  line  motion  along  the  x-axis.  These  results  will  be  confirmed  with 
the  stability  analysis  of  the  next  chapter. 
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Figure  4. Ship  Offset  vs.  Time 


Figure  5.  Leading  Ship  Rudder  Angle  vs.  Time 
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Figure  6.  Towline  Angle  vs.  Time 


Figure  7.  Leading  Ship  Heading  Angle  vs.  Time 
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Figure  8.  Trailing  Ship  Heading  Angle  vs.  Time 

2.  Unstable  Condition 

The  figures  shown  below  were  obtained  for  the  following  parameters: 
Tension=0.005 
Length=1.2 
Time  Constant=T  .2 

The  results  clearly  show  an  oscillatory  divergence,  and  the  system  is 
unstable  despite  the  stabilizing  effect  of  the  leading  ship  control  law.  The 
stability  analysis  of  the  following  chapter  will  explain  these  results. 
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Figure  9.  Ship  Offsets  vs.  Time 


Figure  10.  Leading  Ship  Rudder  Angle  vs.  Time 
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Figure  11.  Towline  Angle  vs.  Time 


Figure  12.  Leading  Ship  Fleading  Angle  vs.  Time 
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Figure  13.  Trailing  Ship  Heading  Angle  vs.  Time 
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III.  STABILITY  ANALYSIS 


A.  LINEARIZATION 

Since  the  equations  of  motion  have  changed  by  the  inclusion  effect  of  rudder,  we 
have  to  again  linearize  the  equations  to  see  whether  the  system  is  stable  or  not.  We  will  use 

the  same  assumptions  that  are  mentioned  before  in  Chapter  2  Section  2. 

When  we  linearize  equations  (2),  (3),  (8),  (9),  (13),  (42),  (43)  we  get 


(m2  ~  Yv2  >2  -  Y,2r2  =  ~m2r2  +  Yv2v2  +  Yr2r2  -T(y/2+y)  (45) 

{I zi  ~  i-2 )c  ~Ni2v2  —  Nv2v2  +  N r2r2  —Txp2(l//2  +  y)  (46) 

K  -  n,  )1>1  -  Yr\r\  =  -"Vi  +  Yv\v\  +  Yrirx  +  T{y/X  +  y)  +  YsS  (47) 

(4i  -  *  «  >i  -  ^ vj  =  Nvlv,  +  Nrlr ,  -  (^  +  y)  +  (48) 


Vj  =  ?/,  +  Vj 


(49) 


y  =v2+v2 


(50) 


sin /  = 


y2  +xP2y/2  -U 
/ 


(51) 


We  can  convert  the  equations  of  motion  into  a  matrix  form  in  order  to  make  it 
easier  study  the  stability  of  the  system.  Since  we  changed  only  the  equations  of  the  first 
ship,  equations  (24),  (25)  are  the  same,  but  equations  (26)  and  (27)  are  different 


V  A2my2  ^2 vrr2  ~^^2v(^2  If )  (24) 

C  —  ^2 rvV2  ^2rrr2  ^2r  ^ 2  7 )  (25) 

Vi  =  4wh  +  +  K  Wx  +y)+ Clvs  (52) 

fi  =  4rvVl  +  Krh  +  (^1  +  j)  +  (53) 
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The  coefficients  are  the  same  as  in  Chapter  2  Section  2,  with  the  following 
exceptions: 


c„ =[(/, -Nn)Ys  +  N,Y„M'", -n,X/a  (54) 

Clr  =  [YsNa  +  Ns (m,  -  K„ )] H-  [(,»,,  -  7,,  )(/„  -  N„ )  -  JV(17„ ]  (55) 


As  mentioned  in  Chapter  1,  all  hydrodynamic  coefficients  are  taken  from 
reference  [2]  with  the  exception  of  Ys  and  Ns .  The  matlab  program  shown  in  Appendix 

F  was  used  to  calculate  these  coefficients.  Since  no  data  were  available  on  the  SLICE 
rudder  design,  we  made  an  assumption  of  a  ship  turning  radius  of  three  ship  lengths  under 
fifteen  degrees  of  rudder,  and  calculated  the  necessary  values  of  the  rudder  hydrodynamic 
coefficients  to  achieve  that  turning  radius. 

As  mentioned  before,  the  rudder  feedback  coefficients  ky,kY,  kr and  ky  are 

calculated  by  standard  pole  placement  techniques.  The  matlab  program  shown  in 
Appendix  G  performs  this  calculation. 

Now  we  can  set  our  matrix  to  analysis  the  stability  of  the  system.  The  new  system 
matrix  is  8x8  and  is  written  as  shown  below 


1)2 1 

V2 

c 

V2 

A 

A 

K 

y2 

=  U] 

rx 

y2 

Ti 

y  1 

¥2 

¥2 

¥i\ 

¥x_ 

where 


(56) 
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0 


M= 


An.  ~  GA 

0 

1 

0 

o 


A, 

o 

o 

0 

1 


-c,A 


A 

/ 

_A 

/ 

- 

B2vXp2 

AA  l 

^2vXp\ 

l 

A 

A 

B  1  Bl'Xpl 

A  rX  p\ 

/ 

/ 

2r  l 

l 

A. 

_  A 

-  -  Clvky 

B\vXp2 

Bu.  +  KXp' 

/ 

/ 

l 

lv  l 

A 

_  A 

r  k 

A  X  p2 

B  +KXp' 

/ 

/ 

^1  rKy 

l 

K+  l 

0 

0 

l 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

-Ci  A 


-A  A 


(57) 


After  fonning  the  matrix,  we  used  the  matlab  code  in  Appendix  H  to  find  the 
eigenvalues  of  the  matrix,  which  will  dictate  whether  the  system  is  stable  or  not.  When 
we  run  the  code,  we  experienced  that  the  stability  of  the  system  changes  according  to 
different  values  of  tension,  length  and  the  time  constant.  A  systematic  series  of  runs  was 
conducted  in  order  to  ascertain  the  effects  of  these  parameters. 


B.  DEGREE  OF  STABILITY 

The  critical  eigenvalue  of  the  8x8  system;  i.e.,  the  eigenvalue  that  dictates 
stability  or  instability  is  the  one  that  has  the  largest  real  part.  We  define  this  as  the  degree 
of  stability  of  the  system.  A  positive  value  indicates  instability  while  a  negative  value 
indicates  stability.  A  systematic  series  of  runs  was  then  conducted  where  the  towline 
length,  the  towline  tension,  and  the  leading  ship  control  time  constant  were  varied.  For 
this  analysis  we  used  the  matlab  code  in  Appendix  I  for  varying  length  of  the  connection 
between  the  ships,  the  code  in  Appendix  J  for  varying  tension  on  the  connection,  and  the 
code  in  Appendix  K  for  varying  control  time  constant. 

The  results  for  different  tow  lengths  are  shown  in  Appendix  A,  for  different 
tensions  in  Appendix  B,  and  for  different  time  constants  in  Appendix  C.  Based  on  these 
results  we  can  draw  the  following  conclusions: 

1.  Relatively  large  time  constants  (a  slow  control  law)  cannot  guarantee  stability 
under  towing.  Sufficiently  fast  control  laws  may,  depending  on  tension  and 
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towline  length,  stabilize  the  complete  system.  Therefore,  towing  stability  must  be 
a  consideration  during  leading  ship  control  law  design. 

2.  Sufficiently  large  values  of  the  tension  (i.e.;  trailing  ship  resistance)  may  stabilize 
the  complete  two-body  system.  This  result  is  consistent  with  earlier  observations 
reported  in  Ref.  [3]. 

3.  For  very  large  values  of  the  towline  length,  the  effect  of  the  leading  ship  control 
appears  to  be  small.  Therefore,  the  previously  reported  results  in  the  literature  [3] 
can  be  considered  as  a  large  towline  length  of  the  two-body  model  developed  in 
this  work. 


C.  REGIONS  OF  STABILITY 

The  previous  results  can  be  summarized  in  a  single  graph  designating  regions  of 
stability  and  instability.  This  graph  is  shown  in  Figure  14  which  was  produced  by 
utilizing  the  Matlab  code  shown  in  Appendix  L.  The  graph  shows  the  critical  value  of  the 
control  time  constant  for  stability  versus  the  towline  length,  and  is  parameterized  by  the 
tension  in  the  towline.  Combinations  of  (Tc,L,T)  below  the  corresponding  curve  will 
yield  stable  response,  while  combinations  that  are  located  above  the  curve  will  result  in 
system  instability.  A  comparison  between  the  values  of  the  parameters  used  in  the 
previous  two  simulation  cases  and  the  results  shown  in  Figure  14,  demonstrates  the 
agreement  between  numerical  integrations  and  the  theoretical  predictions  of  the  behavior 
of  the  system  based  on  our  stability  analysis. 

Finally,  we  present  a  comparison  between  our  results  and  the  results  presented  in 
Ref.  [3],  which  were  based  on  the  4th  order  system.  In  reference  [3],  the  first  necessary 
stability  criterion  based  on  the  trailing  ship  is 


x 


p 


(58) 


where  xp  is  the  non-dimensional  length  between  the  center  of  the  ship  and  the 
towline  connection  point. 


22 


Different  T  Values 


Figure  14.  The  Stability  Region 


The  second  necessary  criterion  is 


T>-°^-  (59) 

ax 

where  T  is  the  tension  of  the  towed  ship.  If  the  above  tension  is  negative  then  the 
criterion  in  equation  (59)  is  inactive  and  is,  therefore,  automatically  satisfied. 

The  parameters  in  the  above  equations  are: 


a2  =  ( B0ClDl  - A0D12)+-(BqC1D2  +  BaC2D]  - 2Af)D]D1 ) 
+ —  [bqC2D2  —  AaD2  ) 

and 

a2  =  B0C0D]  +-(b0C0D2  -B^eJ 


(60) 


(61) 


The  coefficients  used  in  equations  (60),  and  (61)  are  as  follows 
A0=(Yv-m)(Nr-fz)-NvYr  (62) 
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(63) 


B0  =  {V,  - m)Nr  +  F  (A,,  - Iz )+  Nv (m  -Yr)~  YrNv 


YvNr+Nv{m-Yr) 

(64) 

(m-Y.)xp  +  Av 

(65) 

-(A,  -I:)+  (m  -  Yv  )xp2  +  (AT  +  Y. ) 

(66) 

Nv~Yvxp 

(67) 

-Yvxp2  +  Nvxp  +  (: Yr  -  7,  )xp  -  Nr  +  A, 

(68) 

0 

(69) 

Using  the  hydrodynamic  coefficients  of  the  Kaimalino  estimated  in  [2],  and 
presented  in  the  Matlab  code  in  several  of  the  Appendices  in  this  thesis,  we  can  see  that 
x  is  always  greater  than  the  ratio  of  Nv  to  F  .  Furthermore,  the  critical  tension  is  always 

negative  as  is  estimated  by  the  Matlab  code  shown  in  Appendix  M.  As  a  result,  the 
previous  stability  criteria  would  have  predicted  a  system  that  would  always  be  stable  and 
would  have  missed  the  stability  and  instability  regions  shown  in  Figure  14.  Therefore,  in 
towing  stability  studies,  the  designer  should  consider  not  only  the  effect  of  the  trailing 
ship,  but  also  the  effect  of  the  leading  ship. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

This  thesis  presented  a  comprehensive  study  of  the  stability  properties  of  the  two- 
body  ship  towing  problem.  The  results  of  this  study  can  help  in  establishing  rational 
guidelines  which  should  be  followed  in  order  to  ensure  stability  and  safety  of  close 
proximity  ship  towing  operations.  The  main  results  from  this  study  can  be  summarized  as 
follows: 

1.  Relatively  large  time  constants  (a  slow  control  law)  cannot  guarantee  stability 
under  towing.  Sufficiently  fast  control  laws  may,  depending  on  tension  and 
towline  length,  stabilize  the  complete  system.  Therefore,  towing  stability  must  be 
a  consideration  during  leading  ship  control  law  design.  This  is  in  contrast  with 
previous  studies  where  the  dynamics  of  the  leading  ship  were  not  considered 
important  in  the  overall  analysis. 

2.  Sufficiently  large  values  of  the  tension  (i.e.;  trailing  ship  resistance)  may  stabilize 
the  complete  two-body  system.  This  result  is  consistent  with  earlier  observations 
reported  in  the  literature. 

3.  For  very  large  values  of  the  towline  length,  the  effect  of  the  leading  ship  control 
becomes  less  pronounced.  Therefore,  previously  reported  results  in  the  literature 
can  be  considered  as  a  large  towline  length  of  the  two-body  model  developed  in 
this  work. 
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B.  RECOMMENDATIONS 

Recommendations  for  continuing  studies  are  the  following: 

1.  Study  the  effect  of  additional  geometrical  parameters  such  as  attachment  point 
location  on  the  regions  of  stability  and  instability. 

2.  Perform  a  nonlinear  analysis  in  order  to  analyze  the  mechanism  of  the  loss  of 
stability.  Current  results  indicate  that  the  predominant  mechanism  is  through  the 
generation  of  limit  cycles  (periodic  solutions)  but  this  needs  to  be  substantiated 
further. 
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APPENDIX  A.  DEGREE  OF  STABILITY  VS.  TOWLINE  LENGTH 


Constant  T=0.001 


Fig  A.l 

Constant  T=0  003 


Fig  A.2 
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Degree  of  Stability  Degree  of  Stabilit 


Constant  T=0.005 


Fig  A. 3 

Constant  T=0.0075 


Fig  A.4 
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Constant  T=0.01 


Fig  A. 5 


Constant  Tc=0.1 


Fig  A. 6 
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Constant  Tc=0.3 


Fig  A. 7 


Constant  Tc=0  5 


Fig  A. 8 
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Constant  Tc=0.75 


Fig  A. 9 


Constant  Tc=1.0 


Fig  A.  10 
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Constant  Tc=1.5 


Fig  A.  11 


Constant  Tc=2.0 
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Fig  A.  12 
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Degree  of  Stability  Degree  of  Stability 


APPENDIX  B.  DEGREE  OF  STABILITY  VS.  TENSION 


Constant  Tc=0.1 


FigB.l 


Constant  Tc=0.3 


Fig  B.2 
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Constant  Tc=0.5 


Fig  B.3 


Constant  Tc=0.7 


Fig  B.4 
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Constant  Tc=1 .0 


Fig  B.5 


Constant  Tc=1 .25 


Fig  B.6 
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Degree  of  Stability 


Constant  Tc=1 .50 


Fig  B.7 

Constant  Tc=1.75 


Fig  B.8 
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Degree  of  Stability 


Constant  Tc=2.0 


Fig  B.9 


Constant  L=0.1 


Fig  B.  10 
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Constant  L=0.3 


Fig  B.  11 


Constant  L=0.5 


Fig  B.  12 
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Degree  of  Stability  ,  Degree  of  Stability 


Constant  L=0.7 


Fig  B.  13 


Constant  L=1.0 


Fig  B. 14 
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Constant  L=1 .25 


Fig  B. 15 


Constant  L=1 .50 


Fig  B.  16 
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Degree  of  Stability  Degree  of  Stability 


Constant  L— 1 .75 


Fig  B.  17 


Constant  L=2.0 


Fig  B. 18 
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APPENDIX  C.  DEGREE  OF  STABILITY  VS.  TIME  CONSTANT 


Constant  T=0  001 


FigC.l 


Constant  T=0.002 


FigC.2 
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Degree  of  Stability  Degree  of  Stabilit 


Constant  T=0.003 


FigC.3 


Constant  T=0.004 


Fig  C.4 
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Constant  T=0.005 


Fig  C.5 


Constant  T=0  006 


Fig  C.6 
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Constant  T=0  007 


Fig  C.7 


Constant  T=0  008 


Fig  C.8 
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Constant  T=0.009 


Fig  C.9 


Constant  T=0.01 


Fig  C.  10 
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Fig  C.  11 
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0.15  - 


0.05 


0.1  - 


-0.05 


-0.1  - 


-0.15  - 


-0.2 


-  T=0.001 

-  T=0.002 

-  T=0.003 

T=0.004 

-  T=0.005 

-  T=0.006 

-  T =0.007 

T =0.008 
T=0.009 
-  T=0.01 


Fig  C.  12 
48 


0.15 


Constant  L=0.5 


jn 

cc 

w 

o 

CD 

CD 

CD 

CD 

Q 


0.1  - 


-0.15  - 


-0.2 


T =0.001 
T =0.002 
T=0.003 
T=0.004 
T=0.005 
T=0.006 
T =0.007 
T =0.008 
T =0.009 
T=0.01 


0 


0.2  0.4 


0.6  0.8 


1.2  1.4  1.6  1.8 


2 


Fig  C.  13 


0.15 


Constant  L=0.7 


0.05  - 


cs 

(jo 


®  -0.05  - 


-0.1  - 


-0.15 


-0.2 


Fig  C.  14 


49 


Constant  L=1 .0 
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APPENDIX  D.  MATLAB  CODE  FOR  STABILITY  ANALYSIS  WITHOUT 

CONTROL 


%LTJG  Mersin  GOKCE 
%  Thesis  Work 

%  Model  for  coupled  stability  analysis 
%  Initialization 
ml=0. 018078;  m2=0.018; 
ul  =  l ; 
u2  =  l ; 

T=l; 

L=l; 

xpl=0.5;  xp2=0.5; 

lzl=0 .0007;  Iz2=0 . 00069412; 

Yvl=-0 . 07893; 

Yv2=-0 .1183; 

Yrl=-0 . 004044; 

Yr2=-0 . 0042; 

Nvl=-0 .016428; 

Nv2=-0 .0187; 

Nr l=-0 .010332; 

Nr2=-0 . 0176; 

Yvdotl=-0 .051328; 

Yvdot2=-0 .0184; 

Yr dot 1=0 .005617; 

Yrdot2=-0 .0011; 

Nvdotl=-0 .001945; 

Nvdot2=-0 .0008489; 

Nrdotl=-0 . 00564 ; 

Nrdot2=-0 .0090; 

A3=l/L; 

B3=xpl/L; 

C3=xp2/L; 

o, 

o 

A2vv= [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2vr= [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2v= [- ( (Iz2-Nrdot2) *T) - (Yrdot2*T*xp2 ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2rv= [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2rr= [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r= [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ]  ; 

Alvv= [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl-Yrdotl ) ] ; 

Alvr= [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Blv= [ ( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 

Alrv= [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alrr=[ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 
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Blr= [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 

A=zeros ( 7 ) ; 

A ( 1 , 1 ) =A2vv; 

A  ( 1 , 2 ) =A2vr ; 

A (1, 5) =B2v/L; 

A (1,  6) =B2v+ ( (B2v*xp2) /L) ; 

A (1, 7) = (B2v*xpl ) /L; 

A (2, 1) =A2rv; 

A  (2, 2) =A2rr; 

A (2, 5) =B2r/L; 

A(2,6)=B2r+( (B2r*xp2)/L) ; 

A (2, 7) = (B2r*xpl ) /L; 

A (3,  3) =Alvv; 

A (3, 4) =Alvr; 

A (3,  5) =Blv/L; 

A(3,  6)  =  (Blv*xp2) /L; 

A (3,  7) =Blv+ (  (Blv*xpl) /L) ; 

A ( 4 , 3 ) =Alrv; 

A (4, 4) =Alrr; 

A  (4, 5) =Blr/L; 

A  (4, 6)  =  (Blr*xp2) /L; 

A  (4, 7) =Blr+ (  (Blr*xpl) /L) ; 

A (5,  1) =1; 

A (5,  3) =-l; 

A (5,  6) =1; 

A (5,  7) =-l; 

A (6,  2) =1; 

A ( 7 , 4 ) =1 ; 

eig (A) 
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APPENDIX  E.  MATLAB  CODE  TO  PROVE  ZERO  EIGENVALUE 


syms  al  bl  el  fl  gl 
syms  a2  b2  e2  f2  g2 
syms  c3  d3  e3  f3  g3 
syms  c4  d4  e4  f4  g4 
syms  h 

A= [ al  bl  0  0  el  fl  gl; 
a2  b2  0  0  e2  f2  g2; 

0  0  c3  d3  e3  f3  g3; 

0  0  c4  d4  e4  f4  g4; 

10-1001-1; 

0  1  0  0  0  0  0; 

0001000]; 

B=eye  (7)  . *h; 

C=det (A-B) ; 
h=0 ; 

D=subs (C) ; 

syms  ml  m2  ul  u2  T  L 

syms  xpl  xp2  Izl  Iz2 

syms  Yvl  Yv2  Yrl  Yr2 

syms  Nvl  Nv2  Nrl  Nr2 

syms  Yvdotl  Yvdot2  Yrdotl  Yrdot2 

syms  Nvdotl  Nvdot2  Nrdotl  Nrdot2 

syms  A2vv  A2vr  B2v 

syms  A2rv  A2rr  B2r 

syms  Alvv  Alvr  Blv 

syms  Alrv  Alrr  Blr 

L=l; 

xpl=0 . 5 ; 
xp2=0 . 5 ; 

A2vv= [ ( ( Iz2-Nrdot2 ) *Yv2 ) + (Nv2*Yrdot2 ) ] / [ ( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) 
(Nvdot2*Yrdot2 ) ]  ; 

A2vr= [ (Yr2-m2*u2) * (Iz2-Nrdot2) + (Nr2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2v=- [ ( (Iz2-Nrdot2) *T) + (Yrdot2*T*xp2 ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) 
(Nvdot2*Yrdot2 ) ]  ; 

A2rv= [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) - 
(Nvdot2*Yrdot2 ) ] ; 

A2rr= [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2 
Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r=- [ (T*Nvdot2 ) + (T*xp2* (m2-Yvdot2 ))]/[( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) - 
(Nvdot2*Yrdot2 ) ] ; 

Alvv= [ ( (Izl-Nrdotl) *Yvl) + (Nvl* Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 
(Nvdotl-Yrdotl ) ] ; 

Alvr= [ ( Yr l-ml*ul ) * (Izl-Nrdotl ) + (Nrl*Yrdot 1 ) ] / [ ( (ml-Yvdotl ) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Blv= [ ( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alrv= [ (Yvl*Nvdotl ) + (Nvl* (ml-Yvdotl ))]/[( (ml-Yvdotl ) * ( Izl-Nrdotl ) ) - 
(Nvdotl*Yrdotl ) ]  ; 

Alrr=[ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl 
Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Blr= [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 


55 


al=A2vv; 
bl=A2vr ; 
el=B2v/L; 

fl=B2v+ ( (B2v*xp2) /L) ; 

gl= (B2v*xpl) /L; 

a2=A2rv; 

b2=A2rr ; 

e2=B2r/L; 

f2=B2r+ ( (B2r*xp2) /L) ; 

g2= (B2r*xpl) /L; 

c3=Alvv; 

d3=Alvr ; 

e3=Blv/L; 

f3= (Blv*xp2) /L; 

g3=Blv+ ( (Blv*xpl) /L) ; 

c4=Alrv; 

d4=Alrr ; 

e4=Blr/L; 

f4= (Blr*xp2) /L; 

g4=Blr+ ( (Blr*xpl) /L) ; 

E=subs (D) 
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APPENDIX  F.  MATLAB  CODE  TO  FIND  HYDRODYNAMIC  COEFFICIENTS 


%  to  find  Ydel  and  Ndel 
%  Ndel=-0 . 5Ydel 
%  R=1 / (K*del) 

%  K= [ (Nvl*Ydel) - (Yvl*Ndel) ] / [ (Yvl*Nrl) - (Nvl* (Yrl- (ml*ul) ) ) ] 
%  R=3  (Assumption, tactical  diameter) 

%  del=15  degree=  0.262  radian  (Assumption, rudder) 
ml=. 018078; 
ul  =  l ; 

Yvl=-0 . 07893; 

Yrl=-0 . 004044; 

Nvl=-0 . 016428; 

Nr l=-0 .010332; 

R=3; 

del=-0 .262 

Ydel= [Yvl*Nrl- (Nvl* (Yrl-ml) ) ] / [ (Nvl+ (0 . 5*Yvl) ) *R*del] 
Ndel=-0 . 5*Ydel 
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APPENDIX  G.  MATLAB  CODE  TO  FIND  THE  COEFFICIENTS  OF  DELTA 

%  STEP  2 

%  To  Find  Coefficients  of  Delta 
ml=. 018078; 
ul  =  l ; 

T=1 ; 

L=l; 

xpl=0 . 5 ; 

Izl=. 0007; 

Yvl=-0 . 07893; 

Yrl=— 0 . 004044; 

Nvl=-0 . 016428; 

Nr l=-0 .010332; 

Yvdotl=-0 .051328; 

Yr dot 1=0 .005617; 

Nvdotl=-0 .001945; 

Nrdotl=-0 . 00564 ; 

Ydel=0 . 0103; 

Ndel=-0 . 0051 ; 

o, 

o 

"6 

Alvv= [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Alvr= [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / ( ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Blv= [ ( ( Izl-Nrdotl ) *T) - (Yrdotl*T*xpl ) ] / [ ( (ml-Yvdotl ) * ( Izl-Nrdotl ) ) - 
(Nvdotl*Yrdotl ) ] ; 

Clv= [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 

Alrv= [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 

Alrr=[ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 
Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Blr= [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ] ; 

Clr= [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) - 
(Nvdotl*Yrdotl ) ]  ; 

A=zeros ( 4 ) ; 

A ( 1 , 3 ) =1 ; 

A (2, 2) =Alvv; 

A (2, 3) =Alvr; 

A (3,  2) =Alrv; 

A (3,  3) =Alrr; 

A ( 4 , 1 ) =1 ; 

A ( 4 , 2 ) =1 ; 

B=zeros (4,1); 

B (2,  1)  =Clv; 

B (3,  1) =Clr; 

poles= [-0.1  -0.11  -0.12  -0.13]; 
k=place (A, B, poles ) 
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APPENDIX  H.  MATLAB  CODE  FOR  STABILITY  ANALYSIS  AND 

SIMULATION  WITH  CONTROL 

%  Main  program 

%  Eigenvalue  analysis  and  numerical  simulation 


%  Parameters: 
%  T 
%  L 
%  TC 
%  xpl 
amidships ) 

%  xp2 

amidships ) 

%  DeltaT  = 
%  SimTime  = 

o, 

o 


Nondimensional  tension 

Towline  length 

Control  law  time  constant 

=  Attachment  point  (leading  ship  -  positive  forward 

=  Attachment  point  (trailing  ship  -  positive  aft 

Time  step  increment 
Simulation  time 


T 

0 

. 0075; 

L 

0 

.201; 

TC 

0 

•  6; 

bpole  = 

- 

1/TC; 

xpl 

0 

•  5; 

xp2  = 

0 

•  5; 

DeltaT  = 

0 

.  002; 

SimTime= 

2 

o 

o 

NPrint  = 

1 

0; 

o, 

o 

o, 

o 

o, 

o 


Initial  conditions  for  simulation 
yl,y2  must  be  consistent  with  L 


every 


(NPrint ) 


points 


of 

of 


psil_old=0 ; vl_old=0 ; rl_old=0 ; yl_old=0 .01; 
psi2_old=0 ; v2_old=0 ; r2_old=0 ; y2_old=0 .00; 

o, 

o 


%  Constants 


"6 

u  1 

=  i; 

u2 

=  1; 

ml 

=  0.018078; 

m2 

=  0.018; 

Izl 

=  0.0007; 

Iz2 

=  0.00069412; 

Yvl 

=  -0.07893; 

Yv2 

=  -0.1183; 

Yrl 

=  -0.004044; 

Yr2 

=  -0.0042; 

Nvl 

=  -0.016428; 

Nv2 

=  -0.0187; 

Nr  1 

=  -0.010332; 

Nr2 

=  -0.0176; 

Yvdotl 

=  -0.051328; 

Yvdot2 

=  -0.0184; 

Yrdotl 

=  0.005617; 

Yrdot2 

=  -0.0011; 

Nvdotl 

=  -0.001945; 

Nvdot2 

=  -0.0008489; 

Nrdotl 

=  -0.00564; 

Nrdot2 

=  -0.0090; 
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Ydel  =  0.0103; 

Ndel  =  -0.0051; 

A3  =  1/L; 

B3  =  xpl/L; 

C3  =  xp2/L; 

D3  =  -1/L; 

"o 

%  Setup  the  matrix  coefficients 

o_ 

o 

A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) 

(Nvdot2*Yrdot2 ) ]  ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2v  =  [- ( ( Iz2-Nrdot2 ) *T) - (Yrdot2*T*xp2 ) ] / [ ( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) 

(Nvdot2*Yrdot2 ) ]  ; 

A2rv  =  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2-Nrdot2) ) 

(Nvdot2*Yrdot2 ) ] ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2 ) - (T*xp2* (m2-Yvdot2 ))]/[( (m2-Yvdot2 ) * ( Iz2-Nrdot2 ) ) 

(Nvdot2*Yrdot2 ) ]  ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ]  ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Blv  =  [ ( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ]  ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ]  ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ] ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) ) ] / ( ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ]  ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl-Nrdotl) ) 

(Nvdotl*Yrdotl ) ]  ; 

o, 

o 

%  Find  control  gains 

o, 

o 


c 

= 

zeros  ( 4 ) ; 

C  ( 1 , 3 ) 

= 

i; 

C (2, 2) 

= 

Alvv; 

C (2, 3) 

= 

Alvr; 

C (3, 2) 

= 

Alrv; 

C  (3,  3) 

= 

Alrr; 

C  ( 4 , 1 ) 

= 

i; 

C  ( 4 , 2  ) 

= 

l; 

D 

= 

zeros  (4,1); 

D (2, 1) 

= 

Civ; 

D (3, 1) 

= 

Clr; 

poles 

= 

[bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 

k 

= 

place  (C, D, poles ) ; 

Kpsi 

= 

k  ( 1 , 1 )  ; 

Kv 

= 

k  ( 1 , 2 )  ; 

Kr 

= 

k  (1, 3)  ; 

Ky 

= 

k  ( 1 , 4 )  ; 

62 


%  A  matrix 


A 

A  ( 1 , 1 )  = 

A  ( 1 , 2  )  = 

A  ( 1 , 5 )  = 

A  ( 1 ,  6 )  = 

A  ( 1 , 7  )  = 

A  ( 1 , 8  )  = 

A (2, 1)  = 

A (2, 2)  = 

A (2, 5)  = 

A (2, 6)  = 

A (2, 7)  = 

A (2, 8)  = 

A  (3,  3)  = 

A  (3,  4)  = 

A  (3,  5)  = 

A  (3,  6)  = 

A  (3,  7)  = 

A  (3,  8)  = 

A  ( 4 , 3 )  = 

A  ( 4 , 4  )  = 

A  ( 4 , 5 )  = 

A  ( 4 ,  6 )  = 

A  ( 4 , 7  )  = 

A  ( 4 , 8  )  = 

A  (5,  1)  = 

A  (5,  7)  = 

A  (6,  3)  = 

A  (  6,  8  )  = 

A  ( 7 , 2 )  = 

A  ( 8 , 4  )  = 

o, 
o 

%  Find  eigenvalues 

o, 

o 

B=eig (A) 

o_ 

o 

%  Start  simulation  (Euler  fixed  step) 

o, 

o 

NT=SimTime/DeltaT ; 
iPrint=l; iStore=l; 
for  i=l : NT, 

o, 

o 

%  Rudder  angle 
"6 

delta  =  -  (  Kpsi*psil_old  +  Kv*vl_old  +  Kr*rl_old  +  Ky*yl_old  ) ; 

o, 

o 

%  Limit  rudder  angle  between  -0.4  and  +0.4  radians 

o, 

o 

if  delta  >  0.4 
delta=0 . 4 ; 

end 

if  delta  <  -0.4 
delta=-0 . 4 ; 

end 
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zeros  ( 8 ) ; 

A2vv; 

A2vr; 

B2v/L; 

-B2v/L; 

B2v+ (  (B2v*xp2) /L) ; 

(B2v*xpl ) /L; 

A2rv; 

A2rr; 

B2r/L; 

-B2r/L; 

B2r+ ( (B2r*xp2) /L) ; 

(B2r*xpl ) /L; 

Alvv- (Clv*Kv) ; 

Alvr- (Clv*Kr) ; 

Blv/L; 

-Blv/L- (Clv*Ky )  ; 

(Blv*xp2 ) /L; 

Blv+ ( (Blv*xpl) /L) - (Clv*Kpsi)  ; 
Alrv- (Clr*Kv)  ; 

Alrr- (Clr*Kr )  ; 

Blr/L; 

-Blr/L- (Clr*Ky ) ; 

(Blr*xp2 ) /L; 

Blr+ ( (Blr*xpl) /L) - (Clr*Kpsi) ; 

l; 

l; 

l; 

i; 

l; 

l; 


Calculate  angle  gamma 


o, 
o 

o, 
o 

gamma=asin ( (y2_old+xp2*sin (psi2_old) -yl_old+xpl *sin (psil_old) ) /L) ; 

o, 
o 

%  Calculate  xdot=f (x) 

"6 

psildot  =  rl_old; 
vldot  =  Alvv*vl_old 

Blv*sin (gamma+psil_old)  ; 

rldot  =  Alrr*rl_old 

Blr*sin (gamma+psil_old) ; 

yldot  =  ul*sin (psil_old)  +  vl_old*cos (psil_old) ; 
psi2dot  =  r2_old; 

v2dot  =  A2vv*v2_old  +  A2vr*r2_old  +  B2v*sin (gamma+psi2_old) ; 
r2dot  =  A2rr*r2_old  +  A2rv*v2_old  +  B2r*sin (gamma+psi2_old) ; 
y2dot  =  u2*sin (psi2_old)  +  v2_old*cos (psi2_old) ; 

o, 
o 

%  Calculate  new  x  =  old  x  +  Dt  *  xdot 


o, 

o 


psil_new 

=  psil_old 

+ 

DeltaT 

■k 

psildot ; 

vl_new 

=  vl_old 

+ 

DeltaT 

•k 

vldot ; 

rl_new 

=  rl_old 

+ 

DeltaT 

•k 

rldot ; 

yl_new 

=  yl_old 

+ 

DeltaT 

•k 

yldot ; 

psi2_new 

=  psi2_old 

+ 

DeltaT 

•k 

psi2dot; 

v2_new 

=  v2_old 

+ 

DeltaT 

•k 

v2dot ; 

r2_new 

=  r2_old 

+ 

DeltaT 

•k 

r2dot ; 

y2_new 

=  y2_old 

+ 

DeltaT 

•k 

y2dot ; 

o, 

o 


%  Ensure  psi  is  between  2pi  and  -2pi 

o, 

o 

if  psil_new  >  (2*pi) 

psil_new=psil_new-2*pi; 

end 

if  psi2_new  >  (2*pi) 

psi2_new=psi2_new-2*pi ; 

end 

if  psil_new  <  (-2*pi) 

psil_new=psil_new+2*pi; 

end 

if  psi2_new  <  (-2*pi) 

psil_new=psi2_new+2*pi; 

end 

o, 

o 

%  Update  x 

o, 

o 

psil_old=psil_new; vl_old=vl_new; rl_old=rl_new; yl_old=yl_new; 
psi2_old=psi2_new; v2_old=v2_new; r2_old=r2_new; y2_old=y2_new; 

o, 

o 

%  Store  results  every  NPrint  time  steps 

"5 

iPrint=iPrint+l ; 
if  iPrint  >  NPrint 

iStore  =  iStore+1; 

time(iStore)  =  i*DeltaT; 
del  (iStore)  =  delta; 
gam (iStore)  =  gamma; 


+  Alvr*rl_old 
+  Alrv*vl_old 


+  Clv*delta  + 
+  Clr*delta  + 
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psil  (iStore) 

=  psil_new; 

psi2  (iStore) 

=  psi2_new; 

vl  (iStore) 

=  vl_new; 

v2 (iStore) 

=  v2_new; 

rl (iStore) 

=  rl_new; 

r2 (iStore) 

=  r2_new; 

yl  (iStore) 

=  yl_new; 

y2 (iStore) 

=  y2_new; 

iPrint 

=  i; 

end 

end 

"6 

%  Results 

o, 

o 


figure  ( 1 ) 

plot (time, yl , time, y2 ) ; legend (' yl 
figure (2 ) 

plot (time, gam* 5 7 .2958) ; xlabel  (  '  t 
figure ( 3 ) 

plot (time, de 1*57. 2958) ; xlabel ( ' t 
figure  ( 4 ) 

plot (time, psil*57.2958) ; xlabel  (  ' t 
figure  ( 5 ) 

plot (time, psi2*57.2958) ; xlabel  (  ' t 


,  ' y2 ' )  ;  xlabel ( ' y ' ) ; ylabel ( ' y ' ) ; grid 
)  ;  ylabel ( ' \gamma  (deg)  ' ) ; grid 


ylabel  ( ' 

1  \delta 

(deg)  '  ) 

;  grid 

; ylabel 

(  ' \psi_l 

(deg) ’ 

1  ) ; grid 

; ylabel 

(  ' \psi_2 

(deg) ’ 

1  ) ; grid 
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APPENDIX  I.  MATLAB  CODE  TO  FIND  THE  DEGREE  OF  STABILITY  FOR 

DIFFERENT  TOWLINE  LENGTHS 

%  Main  program  for  Length 

%  Eigenvalue  analysis  and  numerical  simulation 


%  Parameters: 

%  T  =  Nondimensional  tension 

%  L  =  Towline  length 

%  TC  =  Control  law  time  constant 

%  xpl  =  Attachment  point  (leading  ship  -  positive  forward 

amidships ) 

%  xp2  =  Attachment  point  (trailing  ship  -  positive  aft 

amidships ) 

o, 

o 

%  Constants 

o, 

o 


u  1 

= 

1; 

u2 

= 

i; 

ml 

= 

0.018078; 

m2 

= 

0.018; 

Izl 

= 

0 . 0007; 

Iz2 

= 

0.00069412; 

Yvl 

= 

-0.07893; 

Yv2 

= 

-0 . 1183; 

Yrl 

= 

-0 . 004044; 

Yr2 

= 

-0.0042; 

Nvl 

= 

-0 . 016428; 

Nv2 

= 

-0 . 0187; 

Nr  1 

= 

-0 . 010332; 

Nr2 

= 

-0 . 0176; 

Yvdotl 

= 

-0.051328; 

Yvdot2 

= 

-0 . 0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0 . 0011; 

Nvdotl 

= 

-0 . 001945; 

Nvdot2 

= 

-0 . 0008489; 

Nrdotl 

= 

-0.00564; 

Nrdot2 

= 

-0 . 0090; 

Ydel 

= 

0 . 0103; 

Ndel 

= 

-0 . 0051; 

T 

= 

\ — 1 
O 

o 

TC 

= 

2.0; 

bpole 

= 

-1/TC; 

xpl 

= 

0.5; 

xp2 

= 

0.5; 

index= 

0; 

for  i 

=0 

.1:0.01:2.0; 

index 

=index+l ; 

L 

=  i; 

A3 

=  1/L; 

B3 

=  xpl / L ; 

C3 

=  xp2/L; 

D3 

=  -1/L; 

%  Setup  the  matrix  coefficients 

o, 

o 


of 

of 
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A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ] ; 

B2 v  =  [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rv  =  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 )  ]  ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Blv  =  [( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl) ]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / ( ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml- 

Yvdotl  ) * ( Izl-Nrdotl ) ) - (Nvdotl*Yrdotl ) ] ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) )]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

o, 

o 

%  Find  control  gains 

o, 

o 

C  =  zeros  ( 4 ) ; 

C ( 1 , 3 )  =  1; 

C  (2, 2)  =  Alvv; 

C (2, 3)  =  Alvr; 

C (3, 2)  =  Alrv; 

C (3, 3)  =  Alrr; 

C ( 4 , 1 )  =  1; 

C ( 4 , 2  )  =  1; 

D  =  zeros  (4,1); 

D  (2, 1)  =  Civ; 

D (3, 1)  =  Clr; 

poles  =  [bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 
k  =  place (C, D, poles ) ; 

Kp  s  i  =  k ( 1 , 1 )  ; 

Kv  =  k  ( 1 , 2  )  ; 

Kr  =  k  (1, 3)  ; 

Ky  =  k  (1,4) ; 

o, 

o 

%  A  matrix 


A  =  zeros  ( 8 ) ; 

A  ( 1 , 1 )  =  A2vv; 

A (1, 2)  =  A2vr; 

A (1, 5)  =  B2v/L ; 

A  (1, 6)  =  -B2v/L; 

A(l,7)  =  B2v+ (  (B2v*xp2 ) /L) ; 
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A  ( 1 , 8  ) 
A (2, 1) 
A (2, 2) 
A (2, 5) 
A (2, 6) 
A (2, 7) 
A (2, 8) 
A  (3,  3) 
A  (3,  4) 
A  (3,  5) 
A  (3,  6) 
A  (3,  7) 
A  (3,  8) 
A  ( 4 , 3 ) 
A  ( 4 , 4  ) 
A  ( 4 , 5 ) 
A  ( 4 ,  6) 
A  ( 4 , 7  ) 
A  ( 4 , 8  ) 
A  (5,  1) 
A  (5,  7) 
A  (6,  3) 
A  (  6,  8  ) 
A  ( 7 , 2  ) 
A  ( 8 , 4  ) 

o, 

o 

%  Find 

o, 

o 


=  (B2v*xpl ) /L; 

=  A2rv; 

=  A2rr; 

=  B2r/L; 

=  -B2r/L; 

=  B2r+ ( (B2r*xp2) / L ) ; 

=  (B2r*xpl ) /L; 

=  Alvv- (Clv*Kv) ; 

=  Alvr- (Clv*Kr )  ; 

=  Blv/L ; 

=  -Blv/L- (Clv*Ky ) ; 

=  (Blv*xp2 ) /L; 

=  Blv+ ( (Blv*xpl) / L ) - (Clv*Kpsi)  ; 
=  Alrv- (Clr*Kv)  ; 

=  Alrr- (Clr*Kr )  ; 

=  Blr/L; 

=  -Blr/L- (Clr*Ky )  ; 

=  (Blr*xp2 ) /L; 

=  Blr+ ( (Blr*xpl) / L ) - (Clr*Kpsi)  ; 

=  i; 

=  1; 

=  1; 

=  i; 

=  i; 

=  i; 

eigenvalues 


B=eig (A) ; 

F=real (B) ; 

H (index) =max (F)  ; 
Lv (index) =L; 

end 
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APPENDIX  J.  MATLAB  CODE  TO  FIND  THE  DEGREE  OF  STABILITY  FOR 

DIFFERENT  TENSIONS 

%  Main  program  for  Tension 

%  Eigenvalue  analysis  and  numerical  simulation 


%  Parameters: 


%  T  =  Nondimensional  tension 

%  L  =  Towline  length 

%  TC  =  Control  law  time  constant 

%  xpl  =  Attachment  point  (leading  ship  -  positive  forward 

amidships ) 

%  xp2  =  Attachment  point  (trailing  ship  -  positive  aft 

amidships ) 

o, 

o 

%  Constants 

o, 

o 


u  1 

= 

1; 

u2 

= 

i; 

ml 

= 

0  . 

018078; 

m2 

= 

0  . 

018; 

Izl 

= 

0  . 

0007; 

Iz2 

= 

o . 

00069412; 

Yvl 

= 

-o . 

07893; 

Yv2 

= 

-o . 

1183; 

Yrl 

= 

-o . 

004044; 

Yr2 

= 

-o . 

0042; 

Nvl 

= 

-o . 

016428; 

Nv2 

= 

-o . 

0187; 

Nr  1 

= 

-o . 

010332; 

Nr2 

= 

-o . 

0176; 

Yvdotl 

= 

-0  . 

051328; 

Yvdot2 

= 

-0  . 

0184; 

Yrdotl 

= 

0  . 

005617; 

Yrdot2 

= 

-0  . 

0011; 

Nvdotl 

= 

-0  . 

001945; 

Nvdot2 

= 

-0  . 

0008489; 

Nrdotl 

= 

-0  . 

00564; 

Nrdot2 

= 

-0  . 

0090; 

Ydel 

= 

0  . 

0103; 

Ndel 

= 

-0  . 

0051; 

L 

= 

2 . 0 

t 

TC 

= 

2 . 0 

t 

bpole 

= 

-1/ 

TC; 

xpl 

= 

0 . 5 

r 

xp2 

= 

0 . 5 

r 

A3 

= 

1/L; 

B3 

= 

xpl/L; 

C3 

= 

xp2/L; 

D3 

= 

-1/L; 

index= 

0; 

for  i 

=0 

o 

o 

1:0.001:0 

index=index+l ; 

T  =  i; 

%  Setup  the  matrix  coefficients 

o, 

o 


of 

of 
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A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ] ; 

B2 v  =  [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rv  =  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 )  ]  ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Blv  =  [( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl) ]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / ( ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml- 

Yvdotl  ) * ( Izl-Nrdotl ) ) - (Nvdotl*Yrdotl ) ] ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) )]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

o, 

o 

%  Find  control  gains 

o, 

o 

C  =  zeros  ( 4 ) ; 

C ( 1 , 3 )  =  1; 

C  (2, 2)  =  Alvv; 

C (2, 3)  =  Alvr; 

C (3, 2)  =  Alrv; 

C (3, 3)  =  Alrr; 

C ( 4 , 1 )  =  1; 

C ( 4 , 2  )  =  1; 

D  =  zeros  (4,1); 

D  (2, 1)  =  Civ; 

D (3, 1)  =  Clr; 

poles  =  [bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 
k  =  place (C, D, poles ) ; 

Kp  s  i  =  k ( 1 , 1 )  ; 

Kv  =  k  ( 1 , 2  )  ; 

Kr  =  k  (1, 3)  ; 

Ky  =  k  (1,4) ; 

o, 

o 

%  A  matrix 


A  =  zeros  ( 8 ) ; 

A  ( 1 , 1 )  =  A2vv; 

A (1, 2)  =  A2vr; 

A (1, 5)  =  B2v/L ; 

A  (1, 6)  =  -B2v/L; 

A(l,7)  =  B2v+ (  (B2v*xp2 ) /L) ; 
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A  ( 1 , 8  ) 
A (2, 1) 
A (2, 2) 
A (2, 5) 
A (2, 6) 
A (2, 7) 
A (2, 8) 
A  (3,  3) 
A  (3,  4) 
A  (3,  5) 
A  (3,  6) 
A  (3,  7) 
A  (3,  8) 
A  ( 4 , 3 ) 
A  ( 4 , 4  ) 
A  ( 4 , 5 ) 
A  ( 4 ,  6) 
A  ( 4 , 7  ) 
A  ( 4 , 8  ) 
A  (5,  1) 
A  (5,  7) 
A  (6,  3) 
A  (  6,  8  ) 
A  ( 7 , 2  ) 
A  ( 8 , 4  ) 

o, 

o 

%  Find 

o, 

o 


=  (B2v*xpl ) /L; 

=  A2rv; 

=  A2rr; 

=  B2r/L; 

=  -B2r/L; 

=  B2r+ ( (B2r*xp2) / L ) ; 

=  (B2r*xpl ) /L; 

=  Alvv- (Clv*Kv) ; 

=  Alvr- (Clv*Kr )  ; 

=  Blv/L ; 

=  -Blv/L- (Clv*Ky ) ; 

=  (Blv*xp2 ) /L; 

=  Blv+ ( (Blv*xpl) / L ) - (Clv*Kpsi)  ; 
=  Alrv- (Clr*Kv)  ; 

=  Alrr- (Clr*Kr )  ; 

=  Blr/L; 

=  -Blr/L- (Clr*Ky )  ; 

=  (Blr*xp2 ) /L; 

=  Blr+ ( (Blr*xpl) / L ) - (Clr*Kpsi)  ; 

=  i; 

=  1; 

=  1; 

=  i; 

=  i; 

=  i; 

eigenvalues 


B=eig (A) ; 

F=real (B) ; 

H (index) =max (F)  ; 
Tv ( index) =T ; 

end 
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APPENDIX  K.  MATLAB  CODE  TO  FIND  THE  DEGREE  OF  STABILITY  FOR 

DIFFERENT  TIME  CONSTANTS 

%  Main  program  Time  Constant 

%  Eigenvalue  analysis  and  numerical  simulation 


%  Parameters: 
%  T 
%  L 
%  TC 
%  xpl 
amidships ) 

%  xp2 

amidships ) 

o, 

o 

%  Constants 

o, 

o 


Nondimensional  tension 

Towline  length 

Control  law  time  constant 

=  Attachment  point  (leading  ship  -  positive  forward 
=  Attachment  point  (trailing  ship  -  positive  aft 


u  1 

= 

1 

r 

u2 

= 

1 

r 

ml 

= 

0 

.  018078; 

m2 

= 

0 

.018; 

Izl 

= 

0 

.  0007; 

Iz2 

= 

0 

. 00069412; 

Yvl 

= 

-0 

. 07893; 

Yv2 

= 

-0 

.1183; 

Yrl 

= 

-0 

. 004044; 

Yr2 

= 

-0 

.0042; 

Nvl 

= 

-0 

.016428; 

Nv2 

= 

-0 

.0187; 

Nr  1 

= 

-0 

.010332; 

Nr2 

= 

-0 

.0176; 

Yvdotl 

= 

-0 

. 051328; 

Yvdot2 

= 

-0 

.0184; 

Yrdotl 

= 

0 

.  005617; 

Yrdot2 

= 

-0 

. 0011; 

Nvdotl 

= 

-0 

.  001945; 

Nvdot2 

= 

-0 

.  0008489; 

Nrdotl 

= 

-0 

. 00564; 

Nrdot2 

= 

-0 

. 0090; 

Ydel 

= 

0 

.0103; 

Ndel 

= 

-0 

. 0051; 

T 

= 

0  . 

\ — 1 
o 

L 

= 

2  . 

0; 

xpl 

= 

0  . 

5; 

xp2 

= 

0  . 

5; 

A3 

= 

1 

/L; 

B3 

= 

xpl/L; 

C3 

= 

xp2/L; 

D3 

= 

-1 

/  L; 

index= 

0; 

for  i 

=0 

.  1 

:0. 01:2.0; 

index=index+l ; 

TC  =i ; 

bpole  =  -1/TC; 

%  Setup  the  matrix  coefficients 

o, 

o 


of 

of 
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A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ] ; 

B2 v  =  [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rv  =  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 )  ]  ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Blv  =  [( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl) ]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / ( ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml- 

Yvdotl  ) * ( Izl-Nrdotl ) ) - (Nvdotl*Yrdotl ) ] ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) )]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

o, 

o 

%  Find  control  gains 

o, 

o 

C  =  zeros  ( 4 ) ; 

C ( 1 , 3 )  =  1; 

C  (2, 2)  =  Alvv; 

C (2, 3)  =  Alvr; 

C (3, 2)  =  Alrv; 

C (3, 3)  =  Alrr; 

C ( 4 , 1 )  =  1; 

C ( 4 , 2  )  =  1; 

D  =  zeros  (4,1); 

D  (2, 1)  =  Civ; 

D (3, 1)  =  Clr; 

poles  =  [bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 
k  =  place (C, D, poles ) ; 

Kp  s  i  =  k ( 1 , 1 )  ; 

Kv  =  k  ( 1 , 2  )  ; 

Kr  =  k  (1, 3)  ; 

Ky  =  k  (1,4) ; 

o, 

o 

%  A  matrix 


A  =  zeros  ( 8 ) ; 

A  ( 1 , 1 )  =  A2vv; 

A (1, 2)  =  A2vr; 

A (1, 5)  =  B2v/L ; 

A  (1, 6)  =  -B2v/L; 

A(l,7)  =  B2v+ (  (B2v*xp2 ) /L) ; 
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A  ( 1 , 8  ) 
A (2, 1) 
A (2, 2) 
A (2, 5) 
A (2, 6) 
A (2, 7) 
A (2, 8) 
A  (3,  3) 
A  (3,  4) 
A  (3,  5) 
A  (3,  6) 
A  (3,  7) 
A  (3,  8) 
A  ( 4 , 3 ) 
A  ( 4 , 4  ) 
A  ( 4 , 5 ) 
A  ( 4 ,  6) 
A  ( 4 , 7  ) 
A  ( 4 , 8  ) 
A  (5,  1) 
A  (5,  7) 
A  (6,  3) 
A  (  6,  8  ) 
A  ( 7 , 2  ) 
A  ( 8 , 4  ) 

o, 

o 

%  Find 

o, 

o 


=  (B2v*xpl ) /L; 

=  A2rv; 

=  A2rr; 

=  B2r/L; 

=  -B2r/L; 

=  B2r+ ( (B2r*xp2) / L ) ; 

=  (B2r*xpl ) /L; 

=  Alvv- (Clv*Kv) ; 

=  Alvr- (Clv*Kr )  ; 

=  Blv/L ; 

=  -Blv/L- (Clv*Ky ) ; 

=  (Blv*xp2 ) /L; 

=  Blv+ ( (Blv*xpl) / L ) - (Clv*Kpsi)  ; 
=  Alrv- (Clr*Kv)  ; 

=  Alrr- (Clr*Kr )  ; 

=  Blr/L; 

=  -Blr/L- (Clr*Ky )  ; 

=  (Blr*xp2 ) /L; 

=  Blr+ ( (Blr*xpl) / L ) - (Clr*Kpsi)  ; 

=  i; 

=  1; 

=  1; 

=  i; 

=  i; 

=  i; 

eigenvalues 


B=eig (A) ; 

F=real (B) ; 

H (index) =max (F)  ; 
TCv (index) =TC; 

end 
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APPENDIX  L.  MATLAB  CODE  TO  FIND  THE  REGION  OF  STABILITY 


%  Stability  graph  -  TC  vs .  L  -  constant  T 

o, 

o 


%  Parameters: 


%  T 
%  L 
%  TC 
%  xpl 
amidships ) 
%  xp2 

amidships ) 


=  Nondimensional  tension 
=  Towline  length 
=  Control  law  time  constant 

=  Attachment  point  (leading  ship  -  positive  forward 

=  Attachment  point  (trailing  ship  -  positive  aft 


o, 

o 

%  Constants 

o, 

o 


u  1 

= 

i; 

u2 

= 

i; 

ml 

= 

0 . 018078; 

m2 

= 

0.018; 

Izl 

= 

0.0007; 

Iz2 

= 

0 . 00069412; 

Yvl 

= 

-0 . 07893; 

Yv2 

= 

-0 . 1183; 

Yrl 

= 

-0 . 004044; 

Yr2 

= 

-0 . 0042; 

Nvl 

= 

-0 . 016428; 

Nv2 

= 

-0 . 0187; 

Nr  1 

= 

-0.010332; 

Nr2 

= 

-0.0176; 

Yvdotl 

= 

-0 . 051328; 

Yvdot2 

= 

-0 . 0184; 

Yrdotl 

= 

0.005617; 

Yrdot2 

= 

-0.0011; 

Nvdotl 

= 

-0.001945; 

Nvdot2 

= 

-0 . 0008489; 

Nrdotl 

= 

-0 . 00564; 

Nrdot2 

= 

-0 . 0090; 

Ydel 

= 

0.0103; 

Ndel 

= 

-0.0051; 

indexL 

= 

0; 

index 

= 

0; 

xpl 

= 

0 .5; 

xp2 

= 

0 .5; 

o_ 

o 

%  Enter  constant 


T 


of 

of 


T  =  input ( ' T  =  '  ) ; 

o, 

o 

%  Start  loop  on  length 

"5 

for  iL  =0.1:0.01:2.0; 

indexL  =  indexL+1; 

indexTC  =  0; 

L  =  iL 


L_v (indexL) =  L; 

A3  =  L/L; 
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B3  =  xpl/L; 

C3  =  xp2/L; 

D3  =  — 1/L; 

O, 

o 

%  Loop  on  TC 

o, 

o 

for  iTC  =0.1:0.01:2.0; 
indexTC=indexTC+l ; 

TC  =iTC; 

bpole  =  -1/TC; 

TC_v (indexTC) =TC; 

"5 

%  Setup  the  matrix  coefficients 

o_ 

o 

A2vv  =  [ ( (Iz2-Nrdot2) *Yv2) + (Nv2*Yrdot2) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2vr  =  [ ( (Yr2- (m2*u2) ) * (Iz2-Nrdot2) ) + (Nr2*Yrdot2) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ] ; 

B2v  =  [-((Iz2-Nrdot2)*T)-(Yrdot2*T*xp2)]/[((m2-Yvdot2)*(Iz2- 

Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rv  =  [ (Yv2*Nvdot2) + (Nv2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

A2rr  =  [ ( (Yr2- (m2*u2) ) *Nvdot2) + (Nr2* (m2-Yvdot2) ) ] / [ ( (m2- 

Yvdot2 ) * ( Iz2-Nrdot2 ) ) - (Nvdot2*Yrdot2 ) ]  ; 

B2r  =  [- (T*Nvdot2) - (T*xp2* (m2-Yvdot2) ) ] / [ ( (m2-Yvdot2) * (Iz2- 

Nrdot2) ) - (Nvdot2*Yrdot2 ) ]  ; 

Alvv  =  [ ( (Izl-Nrdotl) *Yvl) + (Nvl*Yrdotl) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Alvr  =  [ (Yrl-ml*ul) * (Izl-Nrdotl) + (Nrl*Yrdotl) ] / [ ( (ml- 

Yvdotl)  * (Izl-Nrdotl) ) - (Nvdotl*Yrdotl ) ]  ; 

Blv  =  [( (Izl-Nrdotl) *T) - (Yrdotl*T*xpl) ]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Civ  =  [ ( (Izl-Nrdotl) *Ydel) + (Ndel*Yrdotl ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Alrv  =  [ (Yvl*Nvdotl) + (Nvl* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ]  ; 

Alrr  =  [ ( (Yrl- (ml*ul) ) *Nvdotl) + (Nrl* (ml-Yvdotl) ) ] / [ ( (ml- 

Yvdotl  ) * ( Izl-Nrdotl ) ) - (Nvdotl*Yrdotl ) ]  ; 

Blr  =  [ (T*Nvdotl) - (T*xpl* (ml-Yvdotl) )]/[( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

Clr  =  [ (Ydel*Nvdotl ) + (Ndel* (ml-Yvdotl) ) ] / [ ( (ml-Yvdotl) * (Izl- 

Nrdotl)  ) - (Nvdotl*Yrdotl ) ] ; 

o, 

o 

%  Find  control  gains 

o, 

o 


c 

=  zeros ( 4 ) ; 

C  (1, 3) 

=  i; 

c  (2,2) 

=  Alvv; 

C (2, 3) 

=  Alvr; 

C (3, 2) 

=  Alrv; 

C (3, 3) 

=  Alrr; 

C  ( 4 , 1 ) 

=  i; 

C  ( 4 , 2 ) 

=  1; 

D 

=  zeros (4,1); 

D (2, 1) 

=  Civ; 

D (3, 1) 

=  Clr; 

poles 

=  [bpole  bpole-0.05  bpole-0.10  bpole-0.15]; 
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k  =  place (C, D, poles ) ; 

Kp  s  i  =  k  ( 1 ,  1 )  ; 

Kv  =  k  ( 1 , 2  )  ; 

Kr  =  k  (1, 3)  ; 

Ky  =  k  ( 1 , 4 )  ; 

o, 

o 

%  A  matrix 
"6 

A  =  zeros  ( 8 ) ; 

A  ( 1 , 1 )  =  A2vv; 

A (1, 2)  =  A2vr; 

A (1, 5)  =  B2v/L ; 

A  (1, 6)  =  -B2v/L; 

A(l,7)  =  B2v+ ( (B2v*xp2 ) /L) ; 

A  ( 1 , 8 )  =  (B2v*xpl)/L; 

A  (2, 1)  =  A2rv; 

A  (2, 2)  =  A2rr; 

A  (2, 5)  =  B2r/L; 

A  (2, 6)  =  -B2r/L; 

A  (2, 7)  =  B2r+ ( (B2r*xp2) /L)  ; 

A  (2, 8)  =  (B2r*xpl)/L; 

A  (3, 3)  =  Alvv- (Clv*Kv)  ; 

A  (3, 4)  =  Alvr- (Clv*Kr)  ; 

A  (3, 5)  =  Blv/L ; 

A  (3, 6)  =  -Blv/L- (Clv*Ky) ; 

A(3,7)  =  (Blv*xp2)/L; 

A  (3, 8)  =  Blv+ ( (Blv*xpl) /L) - (Clv*Kpsi) ; 

A  ( 4 , 3 )  =  Alrv- (Clr*Kv)  ; 

A  ( 4 , 4 )  =  Alrr- (Clr*Kr )  ; 

A  (4, 5)  =  Blr/L; 

A  (4, 6)  =  -Blr/L- (Clr*Ky)  ; 

A  ( 4 , 7 )  =  (Blr*xp2)/L; 

A  ( 4 , 8 )  =  Blr+  (  (Blr*xpl) /L) - (Clr*Kpsi) ; 

A (5, 1)  =  1; 

A  (5,  7)  =  1; 

A  (6,  3)  =  1; 

A ( 6, 8  )  =  1; 

A  ( 7 , 2  )  =  1; 

A  ( 8 , 4  )  =  1; 

O, 

o 

%  Find  eigenvalues 
"6 

B=eig (A) ; 

F=real (B) ; 

H_new=max  (F) ; 

o, 

o 

%  Detect  if  H  changes  its  sign 

o, 

o 

if  indexTC==l 

H_old=H_new; 

else 

PROD=H_new*H_old; 
if  PROD  >  0 
"6 

%  H  does  not  change  sign  -  keep  going 
"6 

H_old=H_new; 
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else 


%  H  changed  its  sign  -  find  critical  TC  by 

interpolation 

"6 

index=index+l ; 

TC_crit (index) =- (H_old*TC_v (indexTC) - 
H_new*TC_v ( indexTC-1 ) ) / (H_new-H_old) ; 

H_old=H_new; 

end 

end 

end 

end 

o, 

o 

%  Plot  (L,  crit  TC)  curve 

o, 

o 

%plot (L_v,  TC_crit )  ,  xlabel ( ' L ' ) , ylabel ( ' { T_C }_{ \rm 
crit } ' ) , title  (  ' Constant  T ' ) , grid 


linear 
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APPENDIX  M.  MATLAB  CODE  TO  EVALUATE  THE  STABILITY  CRITERIA 

BASED  ON  ONLY  TRAILING  SHIP 

%Check  the  Stability  Criteria  in  the  paper 

%For  the  Second  Ship 

%First  Criteria  xp2  >  Nv2/Yv2 

%Second  Criteria  T  >  - (Alfa2/Alfal) 

m2= .018; 

u2  =  l ; 

T=0 .01; 

L=1 . 8541; 
xp2=0 . 5 ; 

Iz2=. 00069412; 

Yv2=-0 . 1183; 

Yr2=-0 . 0042; 

Nv2=-0 . 0187; 

Nr2=-0 . 0176; 

Yvdot2=-0 . 0184 ; 

Yrdot2=-0 .0011; 

Nvdot2=-0 .0008489; 

Nrdot2=-0 .0090; 

Xp2= (Nv2/Yv2 ) 

A= [ ( (Yvdot2-m2 ) * (Nrdot2-Iz2 ) ) - (Nvdot2*Yrdot2 ) ] ;  %A0 
B= [ ( (Yvdot2-m2) *Nr2) + (Yv2* (Nrdot2-Iz2) ) + (Nvdot2* (m2-Yr2) ) - 
(Yrdot2*Nv2 ) ] ;  %B0 

C= [ (Yv2*Nr2 ) + (Nv2* (m2-Yr2) ) ] ;  %C0 
Cl= [  ( (m2-Yvdot2 ) *xp2 ) +Nvdot2 ]  ; 

C2= [- (Nrdot2-Iz2)  +  ( (m2-Yvdot2) * (xp2A2) )  +  (xp2* (Nvdot2+Yrdot2 ) )  ]  ; 

Dl= [Nv2- (Yv2*xp2 ) ] ; 

D2= [ - ( Yv2* (xp2 A2 ) )  +  (Nv2*xp2 )  +  ( ( Yr2-Yvdot2 ) *xp2 ) -Nr2+Nvdot2 ]  ; 

E1=0  ; 

Alfal= [  (  (B*C1*D1) - (A* (D1A2)  )  )  +  (  (1/L) *  (  (B*C1*D2)  +  (B*C2*D1)- 
(2*A*D1*D2 )  )  )  +  (  (1/  (LA2)  ) *  (  (B*C2*D2) - (A* (D2A2)  )  )  )  ] 

Alfa2= [  (B*C*D1 )  +  ( (1/L) * ( (B*C*D2) -  (BA2*E1) ) )  ] 

Tl=— (Alfa2/Alfal ) 
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